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ABSTRACT 

An algorithm is developed to obtain the grid sen- 
sitivity with respect to design parameters for aero- 
dynamic optimization. The procedure is advocating 
a novel (geometrical) parameterization using spline 
functions such as NURBS (Non-Uniform Rational B- 
Splines) for defining the airfoil geometry. An interac- 
tive algebraic grid generation technique is employed 
to generate C-type grids around airfoils. The grid 
sensitivity of the domain with respect to geometric 
design parameters has been obtained by direct differ- 
entiation of the grid equations. A hybrid approach 
is proposed for more geometrically complex configu- 
rations such as a wing or fuselage. The aerodynamic 
sensitivity coefficients are obtained by direct differ- 
entiation of the compressible two-dimensional thin- 
layer Navier-Stokes equations. An optimization pack- 
age has been introduced into the algorithm in order 
to optimize the airfoil surface. Results demonstrate 
a substantially improved design due to maximized 
lift/drag ratio of the airfoil. 

1. INTRODUCTION 

An essential element in design and optimization 
of aerodynamic surfaces is acquiring the sensitivity 
of aerodynamic surface forces with respect to de- 
sign parameters. 1-3 Several methods concerning the 
derivation of sensitivity equations are currently avail- 
able. Among the most frequently mentioned are 
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Direct Differentiation (DD), Adjoint Variable (AV), 
Symbolic Differentiation (SD), Automatic Differenti- 
ation (AD), and Finite Difference (FD). Each tech- 
nique has its own unique characteristics. The Direct 
Differentiation, adopted in this study, has the advan- 
tage of being exact, due to direct differentiation of 
governing equations with respect to design parame- 
ters. There are two basic components in obtaining 
aerodynamic sensitivity. They are: (1) obtaining the 
sensitivity of the governing equations with respect to 
the state variables, and (2) obtaining the sensitiv- 
ity of the grid with respect to the design parame- 
ters. The sensitivity of the state variables with re- 
spect to the design parameters are described by a set 
of linear-algebraic relation. These systems of equa- 
tions can be solved directly by a LU decomposition 
of the coefficient matrix. This direct inversion proce- 
dure becomes extremely expensive as the problem di- 
mension increases. A hybrid approach of an efficient 
banded matrix solver with influence of off-diagonal 
elements iterated can be implemented to overcome 
this difficulty. 2 

After reviewing relevant literature, it is apparent 
that one aspect of aerodynamic sensitivity analysis, 
namely grid sensitivity, has not been investigated ex- 
tensively. The grid sensitivity algorithms in most of 
these studies are based on structural design models. 
Such models, although sufficient for preliminary or 
conceptional design, are not acceptable for detailed 
design analysis. Careless grid sensitivity evaluations, 
would introduce gradient errors within the sensitivity 
module, therefore, infecting the overall optimization 
process. Development of an efficient and reliable grid 
sensitivity module with special emphasis on aerody- 
namic applications appears essential. 

Among two major classes of grid generation sys- 
tems (Algebraic, Differential), algebraic grid genera- 
tion systems are ideally suited for achieving this ob- 
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jective. The explicit formulation, resulting in a fast 
and suitable grid, enables direct differentiation of grid 
coordinates with respect to design parameters. 4,5 The 
underlying effort here is to avoid the time consuming 
and costly numerical differentiation. In addition, the 
analytical derivatives are exact, a desirable feature 
for sensitivity analysis. An important ingredient of 
grid sensitivity is the surface parameterization. Thef- 
most general parameterization would be to specify 
every grid point on the surface as a design parame- 
ter. This, although convenient, is unacceptable due 
to high computational cost. It is essential to keep 
the number of parameters as low as possible to avoid 
a surge on computational expenses. An analytical 
parameterization, may alleviate that problem but it 
suffers from lack of generality. A compromise would 
be using spline functions such as a Bezier or B-Spline 
function to represent the surface. In this manner, 
most aerodynamically inclined surfaces can be repre- 
sented with only a few control (design) parameters. 



Figure 1: Seven control point representation of a 

generic airfoil 


2. SURFACE MODELING AND GRID 
GENERATION 


Among many ideas proposed for generating any 
arbitrary surface, the approximative techniques of 
using spline functions are gaining a wide range of 
popularity. The most commonly used approxima- 
tive representation is the Non- Uniform Rational B- 
Sphne (NURBS) function. They provide a powerful 
geometric tool for representing both analytic shapes 
(conics, quadrics, surfaces of revolution, etc.) and 
free- form surfaces. 6 The surface is influenced by a set 
of control points and weights where unlike interpo- 
lating schemes the control points might not be at the 
surface itself. By changing the control points and 
corresponding weights, the designer can influence the 
surface with a great degree of flexibility without com- 
promising the accuracy of the design. The relation for 
a NURBS curve is 

n 

X(r) = y^iJ iiP (r)D i i = 0,....,n (1) 

2 = 0 
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where X(r) is the vector valued surface coordinate in 
the r-direction, D z are the control points (forming a 
control polygon), coi are weights, Ni jP (r) are the p - th 
degree B-Spline basis function, and iU p (r) are known 
as the Rational basis functions. 

Figure 1 illustrates a seven control point represen- 
tation of a generic airfoil. The points at the leading 



Figure 2: Critical fuselage cross-sections 


and trailing edges are fixed. Two control points at 
the 0% chord are used to affect the bluntness of the 
section. Similar procedure can be applied to other 
airfoil geometries such as NACA four or five digit se- 
ries. The choice for number of control points and 
their locations are best determined using an inverse 
B-Spline interpolation of the initial data. 6 The al- 
gorithm yields a system of linear equations with a 
positive and banded coefficient matrix. Therefore, it 
can be solved safely using techniques such as Gaus- 
sian elemination without pivoting. The procedure 
can be easily extended to cross-sectional configura- 
tions, when critical cross-sections are defined by sev- 
eral circular conic sections, and the intermediate sur- 
faces have been generated using linear interpolation 
as shown in Fig. 2. Increasing the weights would 
deform the circular segments to other conic segments 
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Jacobian 



Figure 3: Sample C-type grid 


(elliptic, parabolic, etc.) as desired for different flight 
regions. In this manner, the number of design param- 
eters can be kept to a minimum, which is an impor- 
tant factor in reducing the optimization costs. 

The algebraic grid generation system, used in this 
study, is an explicit mathematical expression of a 
physical domain as a function of a computational 
domain. A methodology based on separating the 
boundary definition from the interior definition is es- 
tablished. The interior is then defined as a function 
of information on the boundaries such as position, 
surface derivatives, and an independent variable. An 
example of such formulation with first-order surface 
derivatives is called Two-Boundary Grid Generation 
(TBGG) technique. 7 This, matches both the func- 
tion and its derivative at the boundaries. Figure 3 
illustrating the resultant sample grid for the airfoil 
geometry using this technique. 


3. FLOW ANALYSIS AND SENSITIVITY 
EQUATION 


J _ d(CO 

d(x,y)' 

The residual R can be expressed in generalized curvi- 
linear coordinates (£, r /) as, 

df d(G - G v ) 
dr] 

where F and G are the inviscid and G v is the viscous 
fluxes. 

The equations are solved in their conservation form 
using an upwind cell-centered finite- volume formula- 
tion. A third-order accurate upwind biased inviscid 
flux balance is used in both streamwise and normal 
directions. The finite- volume equivalent of second- 
order accurate central differences is used for viscous 
terms. The resulting discretization represents the 
residual, R(Q), at each cell depending locally on val- 
ues of Q at nine neighboring cells such that 


I^Zj(Q) — — 1) Qzj+lj Qzj— 2; Qzj+2; 

Qi-i ,j i Q*+i ,j i Qz— 2jj Q z + 2j )• (5) 

The discretized governing equations are implicitly 
advanced in time using Euler implicit method which 
is unconditionally stable for all time steps according 
to Fourier stability analysis. An iterative approxi- 
mate factorization (AF) algorithm have been chosen 
to advance the solution in time until 

R(Q*) « 0 (6) 

where Q* are the steady-state values of the field vari- 
ables. The boundary conditions are implicitly imple- 
mented within the governing equations. The airfoils 
surface is considered to be impermeable and adia- 
batic. A standard no-slip boundary condition with 
zero surface velocity has been selected. The pres- 
sure at the surface is evaluated using a zeroth-order 
extrapolation from the interior cells. The density is 
then calculated using the state equation. 


3.1 Analysis 

The two-dimensional thin-layer Navier-Stokes 
equations can be represented as 


3.2 Sensitivity 

For a steady-state solution (i.e., t — »■ oo), Eq.(6) is 
reduces to 


1 <9Q 

J dt 


R(Q) 


Q 



( 2 ) 


Here, R is the residual and J is the transformation 


R(Q*(p)j x(p), p) = o (7) 

where the explicit dependency of R on grid and vec- 
tor of parameters P is evident. The parameters P 
control the grid X as well as the solution Q*. The 

fundamental sensitivity equation containing | \ 
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and described by Taylor et al. 2 is obtained by direct 
differentiation of Eq.(7) as 


"<9R" 
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It is important to notice that Eq.(8) is a set of 
linear algebraic equations , and the matrices [^q-J 

and [~^] are well understood. The vector quantity 
{ rP } * s so ^ u ^ on to Eq.(8) given the sensitivity 
of the grid with respect to the parameters, A 

direct chain rule differentiation of results in 


f SX 

\dP 
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where X# designates the boundary coordinates. The 
vector w represents the boundary sensitivity 
which is directly related to boundary parameteriza- 
tion, discussed previously. It has the importance of 
being one of the dominant factors in calculating the 

sensitivity of surface forces needed for optimization 
r aX i 

process. The matrix is responsible for field 

grid sensitivity with respect to boundary coordinates 
and it is related to the rules which govern the grid 
generation algorithm. For algebraic generation sys- 

r aX i 

terns, the primary components of vp , are the 

L a -A -b j 

interpolation functions which distribute the interior 
grid. 

The sensitivity of the grid with respect to the 
vector of design parameters X£> = {X z -, U', <^z} T 

can be obtained by direct differentiation of the grid 
equations. 5 As a consequence of using algebraic grid 
generation technique in which the boundary grid has 
the dominant effect on the interior grid, the bound- 
ary grid sensitivity coefficient would also be essential 
in influencing the interior grid sensitivity coefficient. 
Therefore, evaluation of the surface grid sensitivity 
coefficients are the most important part of the analy- 
sis and are directly dependent on the surface parame- 
terization. For practical purposes, the grid sensitivity 
and orthogonality at the far-held boundary has been 
ignored. 

The how sensitivity coefficient | j> can now be 
directly obtained using the fundamental sensitivity 
equation, Eq.(8), as 
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provided that grid sensitivity, jq is known. The 

Jacobian matrix, [~^] , can be evaluated by differ- 
entiating the discrete residuals with respect to 

T> . 

four vertices of each cell. The quantity can be 

obtained using a full matrix solver to account for all 
the non-zero contributions outside of central band- 
width. This, although convenient, is not practical for 
Navier-Stokes equations due to large storage require- 
ments. An alternative would be the use of a hybrid 
direct solver with conventional relaxation strategy. 2 


3.3 Optimization 

An objective of a multidisciplinary optimization of 
a vehicle design is to extremize a payoff function com- 
bining dependent parameters from several disciplines. 
Most optimization techniques require the sensitivity 
of the payoff function with respect to free parameters 
of the system. For a hxed grid and solution condi- 
tions, the only free parameters are the surface design 
parameters. Therefore, the sensitivity of the payoff 
function with respect to design parameters is needed. 
The optimization problem is based on the method of 
feasible directions and the generalized reduced gradi- 
ent method. 8 This method has the advantage of pro- 
gressing rapidly to a near-optimum design with only 
gradient information of the objective and constrained 
functions required. The problem can be defined as 
finding the vector of design parameters X£>, which 
will minimize the objective function /(X^ ) subjected 
to constraints 

gj (X D )<0 j = l,m (11) 

and 

X' D <X D < X u d (12) 

where superscripts denote the upper and lower 
bounds for each design parameter. The optimization 
process proceeds iteratively as 

X n D = X ^- 1 + 7 S n (13) 

where n is the iteration number, S n the vector of 
search direction, and j a scalar move parameter. The 
first step is to determine a feasible search direction 
S n , and then perform a one-dimensional search in this 
direction to reduce the objective function as much as 
possible, subjected to the constraints. 

The present optimization strategy is based on max- 
imizing the lift coefficient, Cl, in response to surface 
perturbation, subject to pre-determined design con- 
straints. Upper and lower bounds set for each de- 
sign parameter and the sensitivity derivatives of the 

obiective function, d S L , and the constraint, 9 £ D , 

J ’ dXzU ’ aXz)’ 
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Figure 4: Optimization strategy loop 

are obtained as previously described. 2,3 Throughout 
the analysis, the drag coefficient, Cd , is to be no 
greater than the value of the initial design. The strat- 
egy, illustrated in Fig. 4, requires that the grid and 
grid fljnsitivity derivatives be provided dynamically 
during the automated optimization process. 


5. RESULTS AND DISCUSSION 
5.1 Grid Sensitivity 

The grid sensitivity of a generic airfoil with respect 
to design parameters using the NURBS parameter- 
ization is discussed in this section. The geometry, 
as shown in Fig. 1, has seven pre-specffied control 
points. The control points are numbered counter- 
clockwise, starting and ending with control points (0 
and 6), assigned to the tail of the airfoil. A total of 21 
design parameters (i.e., three design parameters per 
control point) available for optimization purpose. De- 
pending on desired accuracy and degree of freedom 
for optimization, the number of design parameters 
could be reduced for each particular problem. For 
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Figure 5: Grid sensitivity with respect to Y\ 


the present case, such reduction is achieved by con- 
sidering fixed weights and chord-length. Out of the 
remaining four control points with two degrees of free- 
dom for each, control points 1 and 5 have been chosen 
as a case study. The number of design parameters is 
now reduced to four with X£> = {Xi , Yi, X5, Y$} T , 
with initial values specified in Fig.l. The non-zero 
contribution to the surface grid sensitivity coefficients 
of these control points are the basis functions Ri^r) 
and Rs^r). Figure 5 illustrates the field-grid sensi- 
tivity with respect to design parameter Y\ when the 
far-held boundary is placed one chord-length away 
from the surface. The sensitivity gradients are re- 
stricted only to the region influenced by the elected 
control point. This locality feature of the NURBS pa- 
rameterization makes it a desirable tool for complex 
design and optimization when only a local perturba- 
tion of the geometry is warranted. Similar results can 
be obtained for design control point 5 where the sen- 
sitivity gradients are restricted to the lower portion 
of domain. 
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Figure 6 : Mach number contours (a = 0, M ^ =0.7) 



Optimization Cycles 


Figure 7: Optimization cycles history 


5.2 Flow Sensitivity and Optimization 

The second phase of the problem is obtaining the 
flow sensitivity coefficients using the previously ob- 
tained grid sensitivity coefficients. In order to achieve 
this, according to Eq.( 8 ), a converged flow field solu- 
tion about a fixed design point should be obtained. 
Thf Computation is performed on a C-type grid com- 
posed of 141 points in the streamwise direction with 
101 points on the airfoil surface, and 31 points in 
the normal direction. The far-held and outer bound- 
ary were placed about 20 chord-length away from the 
airfoils. It is apparent that such a coarse grid is in- 
adequate for capturing the full physics of the viscous 
how over an airfoils. Therefore, it should be under- 
stood that the main objective here is not to produce 
a highly accurate how held solution rather than to 
demonstrate the feasibility of the approach. 

The two-dimensional, compressible, thin-layer 
Navier-Stokes equations are solved for a free stream 
Mach number of M ^ = 0.7, Reynolds number 

Re oo = 10 6 , and angle of attack a — 0°. The solution 
is implicitly advanced in time using local time step- 
ping as a means of promoting convergence toward the 
steady-state. The residual is reduced by ten orders 
of magnitude. All computations are performed on 
NASA Langley’s Cray-2 mainframe with a computa- 
tion cost of 0.1209xl0 -3 CPU seconds/iteration/grid 
point. Figure 6 demonstrates the Mach number con- 
tours of the converged solution with the lift and drag 
coefficients of Cl — 0.402 and Cd — 0.063. Due 
to surface curvature, the how accelerates along the 
the upper surface to supersonic speeds, terminated 
by a weak shock wave behind which it becomes sub- 
sonic. The sensitivity coefficient, | X obtained 


by previously described iterative strategy . 2 The aver- 
age relative error has been reduced by three orders 
of magnitude. The sensitivities of the aerodynamic 
forces, such as drag and lift coefficients with respect 
to design parameters {Xi, Yi, X 5 , Ys} T , are obtained 
and results are presented in Table 1. An inspection 
of Table 1 indicates the substantial inhuence of pa- 
rameters Y\ and Y 5 on the aerodynamic forces acting 
on the surface. The upper and lower bounds for these 
design parameters are assigned as 

0.2 < Xi < 0.7, -0.1 < Yi < 0.5, 

0.2 < X 5 < 0.7, - 0.1 < Y 5 < 0 . 2 . 

The optimum design is achieved after 17 optimiza- 
tion cycles and a total of 8807 Cray -2 CPU seconds. 
These high computational costs make minimizing the 
number of design parameters in optimization cycle es- 
sential. Table 2 highlights the initial and final values 
of lift and drag coefficients with a 208% improvement 
in their ratio. Table 3 represents the initial and op- 
timum design parameters with parameters Y\ and Y 5 
having the largest change as expected. The history of 
design parameters deformation during the optimiza- 
tion cycles appears in Fig. 7, where the oscillatory 
nature of design perturbations during the early cycles 
are clearly visible. Figure 8 compares the original and 
optimum geometry of the airfoil. 

Several observations should be made at this point. 
First, although control points 1 and 5 demonstrated 
to have substantial influence on the design of the air- 
foil, they are not the only control points affecting 
the design. In fact, control points 2 and 4 near the 
nose might have greater affect due to sensitive na- 
ture of lift and drag forces on this region. The choice 
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Table 2 Comparison of initial and optimized 
performance variables 



Figure 8: Original and optimized airfoils 


Performance 

Initial 

Optimum 

Percent 

Variables 

Design 

Design 

Change 

C'l 

0.402 

0.845 

+ 110.1 

Cd 

0.063 

0.043 

-31.7 

Lift/Drag Ratio 

6.38 

19.65 

+208 


Table 3 Comparison of initial and optimized design 
parameters 


Design 

Parameters 

Initial 

Design 

Optimum 

Design 

Percent 

Change 

w 

0.5 

0.374 

-25.2 

V 

0.2 

0.134 

-33 

x 5 

0.5 

0.414 

-17.2 

V 

0.05 

0.069 

+38 


of control points 1 and 5 here was largely based on 
their camber like behavior. A complete design and 
optimization should include all the relevant control 
points (e.g., control points 1, 2, 4, and 5). For ge- 
ometries with large number of control points, in order 
to contain the computational costs within a reason- 
able range, a criteria for selecting the most influen- 
tial control points for optimization purposes should 
be established. This decision could be based on the 
already known sensitivity coefficients, where control 
points having the largest coefficients could be chosen 
as design parameters. Secondly, the optimum airfoil 
of Fig. 8 is only valid for this particular example and 
design range. As a direct consequence of the non- 
linear nature of governing equations and their sensi- 
tivity coefficients, the validity of this optimum design 
would be restricted to a very small range of the orig- 
inal design parameters. The best estimate for this 
range would be the finite- differ enpi step size used to 
confirm the sensitivity coefficients (i.e., 10 -3 or less). 
All the airfoils with the original control points within 
this range should conform to the optimum design of 
Fig. 8, while keeping the grid and flow conditions 
fixed. 


Table 1 Aerodynamic sensitivity coefficients 


Generic Wing-Section 

Direct Differentiation 

Design Parameters: 

dC L 

dXn 

dC D 

dXn 

w 

-0.297 

— 3.63xl0 -2 

V 

-5.107 

0.549 

x 5 

0.15 

— 2.04xl0 -2 

V 

2.609 

0.287 


6. CONCLUSIONS AND 
RECOMMENDATIONS 

An algorithm is developed to obtain the grid sen- 
sitivity with respect to design parameters for aero- 
dynamic optimization. The algebraic Two-Boundary 
Grid Generation (TBGG) scheme has been directly 
differentiated with respect to design parameters. 
This formulation has the benefits of being exact, ef- 
ficient, and inexpensive. The airfoil is defined ge- 
ometrically using the NURBS approximation of the 
surface. A substantial increase in aerodynamic per- 
formance variables enforces the feasibility of this ap- 
proach for high level design and optimization. 

It is evident that grid sensitivity plays a signifi- 
cant role in the aerodynamic optimization process. 
The algebraic grid generation scheme presented here 
is intended to demonstrate the elements involved in 
obtaining the grid sensitivity from an algebraic grid 
generation system. Each grid generation formulation 
requires considerable analytical differentiation with 
respect to parameters which control the boundaries 
as well as the interior grid. It is implied that aero- 
dynamic surfaces, such as the airfoil considered here, 
should be parameterized in terms of design parame- 
ters. Due to the high cost of aerodynamic optimiza- 
tion process, it is imperative to keep the number of 
design parameters as low as possible. Analytical pa- 
rameterization, although facilitates this notion, has 
the disadvantage of being restricted to simple geome- 
tries. A geometric parameterization such as NURBS, 
with local sensitivity, has been advocated for more 
complex geometries. 
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Future investigations should include the implemen- 
tation of present approach using larger grid dimen- 
sions, adequate to resolve full physics of viscous flow 
analysis. A grid optimization mechanism based on 
grid sensitivity coefficients with respect to grid pa- 
rameters should be included in the overall optimiza- 
tion process. An optimized grid applied to present 
geometry, should increase the quality and conver- 
gence rate of flow analysis within optimization cycles. 
Other directions could be establishing a link between 
geometric design parameters (e.g., control points and 
weights) and basic physical design parameters (e.g., 
camber and thickness). This would provide a con- 
sistent model throughout the analysis which could 
easily be modified for optimization. Also, the effects 
of including all the relevant control points on the de- 
sign cycles should be investigated. Another contribu- 
tion would be the extension of the current algorithm 
to three-dimensional space for complex applications. 
For three-dimensional applications, even a geometric 
parameterization of a complete aerodynamic surface 
can require a large number of parameters for its defi- 
nition. A hybrid approach can be selected when cer- 
tain sections or skeleton parts of a surface are speci- 
fied with NURBS and interpolation formulas are used 
for intermediate surfaces. 
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